home <- here::here()# Load librarieslibrary(tidyverse)library(tidylog)library(RCurl)library(gllvm)library(RColorBrewer)library(ggrepel)library(vegan)library(patchwork)library(RCurl)library(devtools)library(ggpubr)library(janitor)library(brms)library(modelr)library(tidybayes)library(ggstats)library(broom)library(broom.mixed)# Source map-plotsource_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
Reading layer `StatRec_map_Areas_Full_20170124' from data source
`/Users/maxlindmark/Dropbox/Max work/R/cod-interactions/data/shapefiles/ICES-StatRec-mapto-ICES-Areas/StatRec_map_Areas_Full_20170124.shp'
using driver `ESRI Shapefile'
Simple feature collection with 11074 features and 17 fields
Geometry type: POLYGON
Dimension: XY
Bounding box: xmin: -44 ymin: 36 xmax: 69 ymax: 85
Geodetic CRS: WGS 84
#source(paste0(home, "R/functions/map-plot.R"))# Source code for rotating ordination plot (to make it work in ggplot)source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/rotate.R")#source(paste0(home, "R/functions/rotate.R"))
Read data
d <-read_csv(paste0(home, "/data/clean/aggregated_stomach_data.csv"))
Diet similarity using bar plots and ordination
Ordination using gllvm
d2 <- d |> dplyr::select(ends_with("_tot"))d2$pred_id <- d$pred_idd2 <- d2 |>left_join(d |> dplyr::select(pred_length_cm, pred_weight_g, pred_id, species), by ="pred_id")
species pred_length_cm n
Length:35 Min. : 9 Min. : 17.0
Class :character 1st Qu.:18 1st Qu.:107.0
Mode :character Median :27 Median :287.0
Mean :27 Mean :268.4
3rd Qu.:35 3rd Qu.:412.0
Max. :49 Max. :508.0
ggplot(t, aes(pred_length_cm, n, color = species)) +geom_point()
# Fit gllvm model with two latend variables and no covariates, meaning it becomes an ordinationset.seed(999)fit <-gllvm(y = d_mat, family ="tweedie", num.lv =2)
Note that, the tweedie family is implemented using the extended variational approximation method.
summarise: now 195 rows and 3 columns, one group variable remaining (name)
ungroup: no grouping variables
ungroup: no grouping variables
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `max_size = as.numeric(substr(predator_length_grp, 5, 6))`.
Caused by warning:
! NAs introduced by coercion
summarise: now 105 rows and 3 columns, one group variable remaining (name)
ungroup: no grouping variables
ungroup: no grouping variables
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `max_size = as.numeric(substr(predator_length_grp, 5, 6))`.
Caused by warning:
! NAs introduced by coercion
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `max_size = as.numeric(substr(predator_length_grp, 5, 6))`.
Caused by warning:
! NAs introduced by coercion
area_plot <- area_plot |>mutate(name =str_to_sentence(name))bar_diet <-ggplot(data = area_plot, aes(x = max_size, y = prop, fill = name, color = name)) +geom_col(width =4.3) +geom_text(data = n_dat, aes(x = max_size, y =1.08, label = n), inherit.aes =FALSE,size =0, color ="white") +geom_text(data = n_dat, aes(x = max_size, y =1.04, label = n), inherit.aes =FALSE,size =2) +facet_wrap(~predator, scales ="free") +scale_fill_manual(values = pal2, name ="") +scale_color_manual(values = pal2, name ="") +coord_cartesian(expand =0) +scale_x_continuous(breaks =seq(0, 100, 5)) +labs(y ="Proportion", x ="Max. predator size in group [cm]") +theme(legend.key.size =unit(0.2, 'cm'),#legend.text = element_text(face = "italic"),legend.position ="bottom",legend.margin =margin(-0.3, 0, 0, 0, unit ="cm"))
prior class coef group resp
(flat) b
(flat) b mean_biom_sc
(flat) b nameflounder_small_cod
(flat) b nameflounder_small_cod:mean_biom_sc
(flat) b namesmall_cod_large_cod
(flat) b namesmall_cod_large_cod:mean_biom_sc
(flat) b
(flat) b nameflounder_small_cod
(flat) b namesmall_cod_large_cod
(flat) b
(flat) b nameflounder_small_cod
(flat) b namesmall_cod_large_cod
student_t(3, 0, 2.5) Intercept
student_t(3, 0, 2.5) Intercept
logistic(0, 1) Intercept
dpar nlpar lb ub source
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
(vectorized)
phi default
phi (vectorized)
phi (vectorized)
zi default
zi (vectorized)
zi (vectorized)
default
phi default
zi default
stancode(fit)
// generated with brms 2.20.4
functions {
/* zero-inflated beta log-PDF of a single response
* Args:
* y: the response value
* mu: mean parameter of the beta distribution
* phi: precision parameter of the beta distribution
* zi: zero-inflation probability
* Returns:
* a scalar to be added to the log posterior
*/
real zero_inflated_beta_lpdf(real y, real mu, real phi, real zi) {
row_vector[2] shape = [mu * phi, (1 - mu) * phi];
if (y == 0) {
return bernoulli_lpmf(1 | zi);
} else {
return bernoulli_lpmf(0 | zi) +
beta_lpdf(y | shape[1], shape[2]);
}
}
/* zero-inflated beta log-PDF of a single response
* logit parameterization of the zero-inflation part
* Args:
* y: the response value
* mu: mean parameter of the beta distribution
* phi: precision parameter of the beta distribution
* zi: linear predictor for zero-inflation part
* Returns:
* a scalar to be added to the log posterior
*/
real zero_inflated_beta_logit_lpdf(real y, real mu, real phi, real zi) {
row_vector[2] shape = [mu * phi, (1 - mu) * phi];
if (y == 0) {
return bernoulli_logit_lpmf(1 | zi);
} else {
return bernoulli_logit_lpmf(0 | zi) +
beta_lpdf(y | shape[1], shape[2]);
}
}
// zero-inflated beta log-CCDF and log-CDF functions
real zero_inflated_beta_lccdf(real y, real mu, real phi, real zi) {
row_vector[2] shape = [mu * phi, (1 - mu) * phi];
return bernoulli_lpmf(0 | zi) + beta_lccdf(y | shape[1], shape[2]);
}
real zero_inflated_beta_lcdf(real y, real mu, real phi, real zi) {
return log1m_exp(zero_inflated_beta_lccdf(y | mu, phi, zi));
}
}
data {
int<lower=1> N; // total number of observations
vector[N] Y; // response variable
int<lower=1> K; // number of population-level effects
matrix[N, K] X; // population-level design matrix
int<lower=1> Kc; // number of population-level effects after centering
int<lower=1> K_phi; // number of population-level effects
matrix[N, K_phi] X_phi; // population-level design matrix
int<lower=1> Kc_phi; // number of population-level effects after centering
int<lower=1> K_zi; // number of population-level effects
matrix[N, K_zi] X_zi; // population-level design matrix
int<lower=1> Kc_zi; // number of population-level effects after centering
int prior_only; // should the likelihood be ignored?
}
transformed data {
matrix[N, Kc] Xc; // centered version of X without an intercept
vector[Kc] means_X; // column means of X before centering
matrix[N, Kc_phi] Xc_phi; // centered version of X_phi without an intercept
vector[Kc_phi] means_X_phi; // column means of X_phi before centering
matrix[N, Kc_zi] Xc_zi; // centered version of X_zi without an intercept
vector[Kc_zi] means_X_zi; // column means of X_zi before centering
for (i in 2:K) {
means_X[i - 1] = mean(X[, i]);
Xc[, i - 1] = X[, i] - means_X[i - 1];
}
for (i in 2:K_phi) {
means_X_phi[i - 1] = mean(X_phi[, i]);
Xc_phi[, i - 1] = X_phi[, i] - means_X_phi[i - 1];
}
for (i in 2:K_zi) {
means_X_zi[i - 1] = mean(X_zi[, i]);
Xc_zi[, i - 1] = X_zi[, i] - means_X_zi[i - 1];
}
}
parameters {
vector[Kc] b; // regression coefficients
real Intercept; // temporary intercept for centered predictors
vector[Kc_phi] b_phi; // regression coefficients
real Intercept_phi; // temporary intercept for centered predictors
vector[Kc_zi] b_zi; // regression coefficients
real Intercept_zi; // temporary intercept for centered predictors
}
transformed parameters {
real lprior = 0; // prior contributions to the log posterior
lprior += student_t_lpdf(Intercept | 3, 0, 2.5);
lprior += student_t_lpdf(Intercept_phi | 3, 0, 2.5);
lprior += logistic_lpdf(Intercept_zi | 0, 1);
}
model {
// likelihood including constants
if (!prior_only) {
// initialize linear predictor term
vector[N] mu = rep_vector(0.0, N);
// initialize linear predictor term
vector[N] phi = rep_vector(0.0, N);
// initialize linear predictor term
vector[N] zi = rep_vector(0.0, N);
mu += Intercept + Xc * b;
phi += Intercept_phi + Xc_phi * b_phi;
zi += Intercept_zi + Xc_zi * b_zi;
mu = inv_logit(mu);
phi = exp(phi);
for (n in 1:N) {
target += zero_inflated_beta_logit_lpdf(Y[n] | mu[n], phi[n], zi[n]);
}
}
// priors including constants
target += lprior;
}
generated quantities {
// actual population-level intercept
real b_Intercept = Intercept - dot_product(means_X, b);
// actual population-level intercept
real b_phi_Intercept = Intercept_phi - dot_product(means_X_phi, b_phi);
// actual population-level intercept
real b_zi_Intercept = Intercept_zi - dot_product(means_X_zi, b_zi);
}
